# Loading required packages
library(dplyr); library(glmmTMB); library(lme4); library(parameters); library(Matrix); library(methods); library(mgcv); library(nlme); library(numDeriv); library(stargazer); library(TMB)
# Loading the data
load("state_level_data.RData")

#############
# TABLE A.2 #
#############
# Making a table of means and standard deviations for state-level (non-standardized) variables
# Variables to use
cols <- c('logsum',"democrat","on_appropriations","distance_from_dw_median","logpop","female","previous_vote_share","seniority","party_leader",
        "freshman", "POPPCT_URBAN", "state.median.household.income")

# Table
stargazer(as.data.frame(senators_appropriations_state[,cols]), 
          covariate.labels=c("Log(State Appropriations Requests + 1)","Democrat (Majority Party Member)", "Member of Appropriations Committee",
                             "Distance from Floor Median","Log State Population", "Senator is a Woman","Previous Vote Share", "Seniority", "Party Leader", 
                             "Freshman Senator", "State Percent Urban Population", "State Median Household Income"),
          title="Summary Statistics of Dataset Analyzing State-Level Outcomes", digits=2, summary.stat= c("mean", "sd", "min","max"),
          notes='\\parbox[t]{1\\textwidth}{\\footnotesize \\textit{Note}: Table presents summary statistics for the (non-standardized) variables used in our state-level analysis.}')

#############
# TABLE C.1 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# FY 2022
# Running the zero inflated model
zim.fy2022 <- glmmTMB(logsum ~ logpop_scaled + on_appropriations +  female + democrat +
                   previous_vote_share_scaled + seniority_scaled + party_leader + distance_from_dw_median_scaled +
                   freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                 zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state[senators_appropriations_state$year==2021,], family = gaussian)
# Model summary and standard errors
summary_zim_fy2022 <- summary(zim.fy2022)
ses_fy2022 <- standard_error(zim.fy2022)

# Save coefficients from the model
cond_fy2022 <- summary_zim_fy2022$coefficients$cond
zi_fy2022 <- summary_zim_fy2022$coefficients$zi
coefs_fy2022 <- c(cond_fy2022[,1])
coefs2_fy2022 <- c(zi_fy2022[,1])

# Save standard errors from the model
ses_cond_fy2022 <- c(ses_fy2022$SE[ses_fy2022$Component=="conditional"])
names(ses_cond_fy2022) <- ses_fy2022$Parameter[ses_fy2022$Component=="conditional"]
ses_zi_fy2022 <- c(ses_fy2022$SE[ses_fy2022$Component=="zero_inflated"])
names(ses_zi_fy2022) <- ses_fy2022$Parameter[ses_fy2022$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_fy2022 <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled +
                  female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                  freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                data = senators_appropriations_state[senators_appropriations_state$year==2021,])

basic_reg2_fy2022 <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                   democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state[senators_appropriations_state$year==2021 & complete.cases(senators_appropriations_state),])

# FY 2023
# Running the zero inflated model
zim.fy2023 <- glmmTMB(logsum ~ logpop_scaled+ on_appropriations +  female + democrat +
                        previous_vote_share_scaled + seniority_scaled + party_leader + distance_from_dw_median_scaled +
                        freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                      zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                      data = senators_appropriations_state[senators_appropriations_state$year==2022,], family = gaussian)
# Model summary and standard errors
summary_zim_fy2023 <- summary(zim.fy2023)
ses_fy2023 <- standard_error(zim.fy2023)

# Save coefficients from the model
cond_fy2023 <- summary_zim_fy2023$coefficients$cond
zi_fy2023 <- summary_zim_fy2023$coefficients$zi
coefs_fy2023 <- c(cond_fy2023[,1])
coefs2_fy2023 <- c(zi_fy2023[,1])

# Save standard errors from the model
ses_cond_fy2023 <- c(ses_fy2023$SE[ses_fy2023$Component=="conditional"])
names(ses_cond_fy2023) <- ses_fy2023$Parameter[ses_fy2023$Component=="conditional"]
ses_zi_fy2023 <- c(ses_fy2023$SE[ses_fy2023$Component=="zero_inflated"])
names(ses_zi_fy2023) <- ses_fy2023$Parameter[ses_fy2023$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_fy2023 <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled  +
                       female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                       freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                     data = senators_appropriations_state[senators_appropriations_state$year==2022,])

basic_reg2_fy2023 <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                        democrat*distance_from_dw_median_scaled,
                      data = senators_appropriations_state[senators_appropriations_state$year==2022 & complete.cases(senators_appropriations_state),])

# Stargazer table of regression results
# coef and se take the coefficients and standard errors from the zero-inflated model
c1.model.list <- list(basic_reg2_fy2022, basic_reg_fy2022, basic_reg2_fy2023, basic_reg_fy2023)
stargazer(c1.model.list,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage 2022","Second Stage 2022","First Stage 2023","Second Stage 2023"),
          dep.var.labels = "Log(State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Median Household Income",
                              "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. All non-binary variables are standardized for ease of comparison. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_state_sepyears",
          digits=3,
          coef=list(coefs2_fy2022, coefs_fy2022, coefs2_fy2023, coefs_fy2023), 
          se=list(ses_zi_fy2022, ses_cond_fy2022, ses_zi_fy2023, ses_cond_fy2023),
          digits.extra = 0,
          title="Predictors of Spending Requests at the State Level Subset by Fiscal Year",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

#############
# TABLE C.4 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# Pooling both years together
senators_appropriations_state_oneyear <- senators_appropriations_state %>%
  group_by(senator,state, logpop, on_appropriations, female, democrat,
           freshman, seniority, party_leader, previous_vote_share, distance_from_dw_median,
           POPPCT_URBAN, state.median.household.income)%>%
  dplyr::summarize(logsum = log(sum(appropriation_sum)+1))

# Re-standardize the non-binary variables (now 100 observations instead of 200)
# Variables to standardize
senators_appropriations_state_scaled <- senators_appropriations_state_oneyear[,c("logpop", "previous_vote_share", "seniority", "distance_from_dw_median", "POPPCT_URBAN", "state.median.household.income")]
# Standardize
senators_appropriations_state_scaled <- as.data.frame(scale(senators_appropriations_state_scaled))
# Rename standardized variables
colnames(senators_appropriations_state_scaled) <- c("logpop_scaled", "previous_vote_share_scaled", "seniority_scaled", "distance_from_dw_median_scaled", "POPPCT_URBAN_scaled", "state.median.household.income_scaled")
# Add the standardized variables to the dataset
senators_appropriations_state_oneyear <- bind_cols(senators_appropriations_state_oneyear, senators_appropriations_state_scaled)

# Pooled
# Running the zero inflated model
zim_pooled <- glmmTMB(logsum ~ logpop_scaled  + on_appropriations +  female + democrat +
                   previous_vote_share_scaled + seniority_scaled + party_leader + distance_from_dw_median_scaled +
                   freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                 zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state_oneyear, family = gaussian)
# Model summary and standard errors
summary_zim_pooled <- summary(zim_pooled)
ses_pooled  <-standard_error(zim_pooled)

# Save coefficients from the model
cond_pooled <- summary_zim_pooled$coefficients$cond
zi_pooled <- summary_zim_pooled$coefficients$zi
coefs_pooled <- c(cond_pooled[,1])
coefs2_pooled <- c(zi_pooled[,1])

# Save standard errors from the model
ses_cond_pooled <- c(ses_pooled$SE[ses_pooled$Component=="conditional"])
names(ses_cond_pooled) <- ses_pooled$Parameter[ses_pooled$Component=="conditional"]
ses_zi_pooled <- c(ses_pooled$SE[ses_pooled$Component=="zero_inflated"])
names(ses_zi_pooled) <- ses_pooled$Parameter[ses_pooled$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_pooled <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled  +
                  female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                  freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled,
                data = senators_appropriations_state_oneyear)

basic_reg2_pooled <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                   democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state_oneyear[complete.cases(senators_appropriations_state_oneyear),])

# Stargazer table of regression results
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_pooled, basic_reg_pooled,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Median Household Income"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. All non-binary variables are standardized for ease of comparison. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$.}",
          label="tab1_state_pooled",
          digits=3,
          coef=list(coefs2_pooled, coefs_pooled), 
          se=list(ses_zi_pooled, ses_cond_pooled),
          digits.extra = 0,
          title="Predictors of Spending Requests at the State Level Pooling Fiscal Years",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

#############
# TABLE C.6 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# Running the zero inflated model
zim_main <- glmmTMB(logsum ~ logpop_scaled  + on_appropriations +  female + democrat + previous_vote_share_scaled + seniority_scaled 
                  + party_leader + distance_from_dw_median_scaled + freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                  zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                  data = senators_appropriations_state, family = gaussian)
# Model summary and standard errors
summary_zim_main <- summary(zim_main)
ses_main <- standard_error(zim_main)

# Save coefficients from the model
cond_main <- summary_zim_main$coefficients$cond
zi_main <- summary_zim_main$coefficients$zi
coefs_main <- c(cond_main[,1])
coefs2_main <- c(zi_main[,1])

# Save standard errors from the model
ses_cond_main <- c(ses_main$SE[ses_main$Component=="conditional"])
names(ses_cond_main) <- ses_main$Parameter[ses_main$Component=="conditional"]
ses_zi_main <- c(ses_main$SE[ses_main$Component=="zero_inflated"])
names(ses_zi_main) <- ses_main$Parameter[ses_main$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_main <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled +
                  female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                  freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                data = senators_appropriations_state)

basic_reg2_main <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                   democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state[complete.cases(senators_appropriations_state),])

# Stargazer table of regression results
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_main, basic_reg_main,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Median Household Income",
                              "Fiscal Year 2023"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. All non-binary variables are standardized for ease of comparison. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_state",
          digits=3,
          coef=list(coefs2_main, coefs_main), 
          se=list(ses_zi_main, ses_cond_main),
          digits.extra = 0,
          title="Predictors of Spending Requests at the State Level",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

#############
# TABLE C.9 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# Set the seed to run the bootstrap and create empty lists to store results
set.seed(02138)
# 20000 bootstraps
nboot <- 20000
results_cond_coefs  <- list()
results_cond_ses  <- list()
results_zi_coefs  <- list()
results_zi_ses  <- list()

# Creating a unique senator variable
senators_appropriations_state$unique_senator_name <- paste0(senators_appropriations_state$senator, senators_appropriations_state$state)

# Bootstrapping at the senator level
ids <- unique(senators_appropriations_state$unique_senator_name)

# Running the bootstrap
# Ensure R and package versions match readme
for (k in 1:nboot) {
  # Sample senators to be used in this bootstrap iteration
  senator_sample <- sample(ids,  length(ids),replace=TRUE)
  
  # All of the senator profiles for this sample
  dataset <- inner_join(tibble(unique_senator_name = senator_sample), senators_appropriations_state, by = "unique_senator_name")
  
  # Then, run the zero-inflated model with this sample
  zim_2 <- glmmTMB(logsum ~ logpop_scaled + on_appropriations +  female + democrat +
                     previous_vote_share_scaled + seniority_scaled + party_leader +
                     distance_from_dw_median_scaled +
                     freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                   zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                   data = dataset, family = gaussian)
  
  summary_zim <- summary(zim_2)
  
  # Storing the results
  results_cond_coefs[[k]]  <- summary_zim$coefficients$cond[,1]
  results_cond_ses[[k]]  <- summary_zim$coefficients$cond[,2]
  results_zi_coefs[[k]]  <- summary_zim$coefficients$zi[,1]
  results_zi_ses[[k]]  <- summary_zim$coefficients$zi[,2]
  
  # Printing progress
  if (k %% 100 == 0) cat(paste(k, "out of", nboot, "samples stored..\n"))
}

# Turn the bootstrap results into our CI estimates
# Omitting models with NA SEs
to.keep <- which(apply(is.na(do.call(rbind.data.frame,results_cond_ses)),1,sum) == 0)

# SD of coefficients to create standard errors for bootstrap
cond_vec_se_alt <- apply(do.call(rbind.data.frame,results_cond_coefs)[to.keep,], 2, sd)
zi_vec_se_alt <- apply(do.call(rbind.data.frame,results_zi_coefs)[to.keep,], 2, sd)

# Variable names for our bootstrapped SEs
names(cond_vec_se_alt) <- c("intercept","logpop_scaled", "on_appropriations", "female", "democrat", "previous_vote_share_scaled", "seniority_scaled", 
                            "party_leader","distance_from_dw_median_scaled", "freshman", "POPPCT_URBAN_scaled", 
                            "state.median.household.income_scaled", "factor(year)2022")
names(zi_vec_se_alt) <- c("intercept","democrat", "on_appropriations", "distance_from_dw_median_scaled", "democrat:distance_from_dw_median_scaled")

# Rounding
cond_vec_se_alt <- round(cond_vec_se_alt,3)
zi_vec_se_alt <- round(zi_vec_se_alt,3)

# Run the zero-inflated model outside of the loop to get coefficients
zim_bootstrap <- glmmTMB(logsum ~ logpop_scaled + on_appropriations +  female + democrat +
                   previous_vote_share_scaled + seniority_scaled + party_leader + distance_from_dw_median_scaled +
                   freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                 zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state, family = gaussian)
# Model summary
summary_zim_bootstrap <- summary(zim_bootstrap)

# Save coefficients from the model
cond_bootstrap <- summary_zim_bootstrap$coefficients$cond
zi_bootstrap <- summary_zim_bootstrap$coefficients$zi
coefs_bootstrap <- c(cond_bootstrap[,1])
coefs2_bootstrap <- c(zi_bootstrap[,1])

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model and the bootstrap above
basic_reg_bootstrap <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled +
                  female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                  freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                data = senators_appropriations_state)

basic_reg2_bootstrap <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                   democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state[complete.cases(senators_appropriations_state),])

# Stargazer table of regression results
# se takes the bootstrapped SEs
stargazer(basic_reg2_bootstrap, basic_reg_bootstrap,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from DW-NOMINATE Median",
                              "Distance from DW-NOMINATE Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Median Household Income",
                              "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. Standard errors come from 19116 bootstraps conducted at the senator level. 
                        We are forced to drop 884 of our initial 20000 bootstrapped iterations due to insufficient coverage for some of our explanatory variables. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_state_boot",
          digits=3,
          coef=list(coefs2_bootstrap, coefs_bootstrap), 
          se=list(zi_vec_se_alt, cond_vec_se_alt),
          digits.extra = 0,
          title="Predictors of Spending Requests at the State Level with Bootstrapped Standard Errors",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.12 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# Running the per capita model
zim.logpercap <- glmmTMB(log_percap ~  logpop_scaled  + on_appropriations +  female + democrat + previous_vote_share_scaled + seniority_scaled 
                         + party_leader + distance_from_dw_median_scaled + freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year),
                         zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                         data = senators_appropriations_state, family = gaussian)
# Model summary and standard errors
summary_zim_logpercap <- summary(zim.logpercap)
ses_logpercap <-standard_error(zim.logpercap)

# Save coefficients from the model
cond_logpercap <- summary_zim_logpercap$coefficients$cond
zi_logpercap <- summary_zim_logpercap$coefficients$zi
coefs_logpercap <- c(cond_logpercap[,1])
coefs2_logpercap <- c(zi_logpercap[,1])

# Save standard errors from the model
ses_cond_logpercap <- c(ses_logpercap$SE[ses_logpercap$Component=="conditional"])
names(ses_cond_logpercap) <- ses_logpercap$Parameter[ses_logpercap$Component=="conditional"]
ses_zi_logpercap <- c(ses_logpercap$SE[ses_logpercap$Component=="zero_inflated"])
names(ses_zi_logpercap) <- ses_logpercap$Parameter[ses_logpercap$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_percap <- lm(log_percap ~  democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled  +
                         female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                         freshman + POPPCT_URBAN_scaled + state.median.household.income_scaled + factor(year), data = senators_appropriations_state)

basic_reg2_percap <- lm(log_percap ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled, 
                        data = senators_appropriations_state[complete.cases(senators_appropriations_state),])

# Stargazer table of regression results
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_percap, basic_reg_percap,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(Per Capita State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Median Household Income",
                              "Fiscal Year 2023"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. All non-binary variables are standardized for ease of comparison. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_statepercap",
          digits=3,
          coef=list(coefs2_logpercap, coefs_logpercap), 
          se=list(ses_zi_logpercap, ses_cond_logpercap),
          digits.extra = 0,
          title="Predictors of Per Capita Spending Requests at the State Level",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.21 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("state_level_data.RData")

# Running the below poverty line model
zim_below_poverty <- glmmTMB(logsum ~ logpop_scaled  + on_appropriations +  female + democrat +
                   previous_vote_share_scaled + seniority_scaled + party_leader + distance_from_dw_median_scaled +
                   freshman + POPPCT_URBAN_scaled + state.pct.poverty_scaled + factor(year),
                 zi = ~  democrat + on_appropriations + distance_from_dw_median_scaled + democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state, family = gaussian)
# Model summary and standard errors
summary_zim_below_poverty <- summary(zim_below_poverty)
ses_below_poverty <-standard_error(zim_below_poverty)
# Save coefficients from the model
cond_below_poverty <- summary_zim_below_poverty$coefficients$cond
zi_below_poverty <- summary_zim_below_poverty$coefficients$zi
coefs_below_poverty <- c(cond_below_poverty[,1])
coefs2_below_poverty <- c(zi_below_poverty[,1])

# Save standard errors from the model
ses_cond_below_poverty <- c(ses_below_poverty$SE[ses_below_poverty$Component=="conditional"])
names(ses_cond_below_poverty) <- ses_below_poverty$Parameter[ses_below_poverty$Component=="conditional"]
ses_zi_below_poverty <- c(ses_below_poverty$SE[ses_below_poverty$Component=="zero_inflated"])
names(ses_zi_below_poverty) <- ses_below_poverty$Parameter[ses_below_poverty$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_below_poverty <- lm(logsum ~ democrat + on_appropriations + distance_from_dw_median_scaled + logpop_scaled  +
                  female + previous_vote_share_scaled + seniority_scaled + party_leader  +
                  freshman + POPPCT_URBAN_scaled + state.pct.poverty_scaled + factor(year),
                data = senators_appropriations_state)

basic_reg2_below_poverty <- lm(logsum ~  democrat + on_appropriations + distance_from_dw_median_scaled +
                   democrat*distance_from_dw_median_scaled,
                 data = senators_appropriations_state[complete.cases(senators_appropriations_state),])

# Stargazer table of regression results
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_below_poverty, basic_reg_below_poverty,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(State Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log State Population",
                              "Senator is a Woman","Previous Vote Share" , "Seniority", "Party Leader", 
                              "Freshman Senator", "State Percent Urban Population", "State Percent Below the Poverty Line",
                              "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's state-level appropriation request behavior. All non-binary variables are standardized for ease of comparison. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_state_poverty",
          digits=3,
          coef=list(coefs2_below_poverty, coefs_below_poverty), 
          se=list(ses_zi_below_poverty, ses_cond_below_poverty),
          digits.extra = 0,
          title="Predictors of Spending Requests at the State Level Replacing Median Income with Percent Below the Poverty Line",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))
